Dynamics of thyroid diseases and thyroid‐axis gland masses

Abstract Thyroid disorders are common and often require lifelong hormone replacement. Treating thyroid disorders involves a fascinating and troublesome delay, in which it takes many weeks for serum thyroid‐stimulating hormone (TSH) concentration to normalize after thyroid hormones return to normal. This delay challenges attempts to stabilize thyroid hormones in millions of patients. Despite its importance, the physiological mechanism for the delay is unclear. Here, we present data on hormone delays from Israeli medical records spanning 46 million life‐years and develop a mathematical model for dynamic compensation in the thyroid axis, which explains the delays. The delays are due to a feedback mechanism in which peripheral thyroid hormones and TSH control the growth of the thyroid and pituitary glands; enlarged or atrophied glands take many weeks to recover upon treatment due to the slow turnover of the tissues. The model explains why thyroid disorders such as Hashimoto's thyroiditis and Graves' disease have both subclinical and clinical states and explains the complex inverse relation between TSH and thyroid hormones. The present model may guide approaches to dynamically adjust the treatment of thyroid disorders.


Appendix Supplementary Text: Mathematical Model for glandmass dynamics in the thyroid axis Model definition
To Model HPT dynamics, we use the following equations: Where the variables are:

Pituitary thyrotroph mass
We summarize the model parameters in the following table: Carrying capacity terms for the thyroid/pituitary gland respectively. Note that when = 0, = 0, this mean that the glands do not have carrying capacities and can grow indefinitely. Otherwise, T,P maximal masses are 1/ , 1/ respectively. TSH-receptor stimulating antibodies 30 External thyroid hormone supply. Michaelis-Menten coefficient for the response function of TH Equation (1) describes TRH ( 1 ) dynamics. TRH production is stimulated by hypothalamic input u which represents the integrated effect of cues including temperature, illness and nutritional states, and is inhibited by thyroid hormones T3 and T4 ( 3 ).
Equation (3) describes thyroid hormone dynamics ( 3 ), which is secreted by the thyroid gland, whose functional mass is T, when stimulated by TSH ( 2 ). 30 represents external thyroid hormone supply. We use this term when simulating levothyroxine treatment of Hashimoto's thyroiditis. Under normal conditions 30 = 0. Ab represents TSH-receptor stimulating antibodies. We use this term when modeling Grave's disease. Under normal conditions Ab = 0. 2 is a Michaelis-Menten coefficient for the response function of TH. We use this parameter to explore the effect of using a MM response function for TH, which is more realistic than a linear function. We found that this choice does not affect the qualitative results of the model. 1 , 2 , 3 are removed at rates 1 , 2 , 3 respectively. Equation (4) for the thyroid functional mass has a removal/turnover termand a mass growth term activated by TSH (x2). (1 − ) is a carrying capacity term, that allows positive growth rate term only for < 1/ . The thyroid growth rate is also affected by activating antibodies Ab if they exist, and is generally described by a MM function with MM coefficient 2 . Equation (5) describes the dynamics of pituitary thyrotroph mass P, which follows a similar balance of growth and removal, where we assume that the main control of mass growth is a suppressive effect by TH ( 3 ).
is the removal rate of thyrotrophs, and the carrying capacity (maximal mass) of thyrotrophs is 1/ .
Under normal conditions, both glands are far from their carrying capacities, ≅ ≅ 0. When modeling Hashimoto's thyroiditis, the relevant carrying capacity is that of the pituitary gland (because the thyroid gland is destructed and is thus small), therefore we approximate ≅ 0. When modeling Graves' disease, the relevant carrying capacity is that of the thyroid gland (the pituitary is small due to the negative TH feedback), therefore we approximate ≅ 0.

Steady state of the simple model
In general, the steady-state of the model cannot be computed explicitly. However, we can compute it for a simple model that represents the healthy state, in which there are no TSHR antibodies and no external thyroid hormone supply ( = 0, 30 = 0), the glands are far from their carrying capacities (i.e. infinite carrying capacity = = 0), and 3 depends approximately linearly on 2 ( 2 = 0). In this case the steady-state is: Therefore, the model equations transform to the dimensionless form: Note that with the choice = 2 = 30 = = = 0 we recover the simple model.
The model shows exact adaptation only when glands are far from their carrying capacities Note that in this simple model, there is exact adaptation for x2 and for x3: the values of these variables depend only on the production and removal rates of the glands. This happens since the glands P and T function as integral feedback controllers, as can be seen in the equations for the simple model below: The gland masses change so that they compensate for any change in the model parameters, except for the production and removal rates of the gland, guaranteeing that 2 = 3 = 1. Note that 2 and 3 were rescaled, so that in terms of the original variables 2 = and 3 = .
Therefore, this simple model cannot account for cases of hypo-or hyperthyroidism.
Adding carrying capacity terms break this perfect adaptation, since the glands can no longer fully compensate for change in the hormone secretion or removal parameters ( Fig EV1). In this case, the glands function as non-ideal integral controllers, where the non-ideal aspect comes naturally from the finite carrying capacity.
Adding a carrying capacity term for the pituitary gland breaks the adaptation of x2 (TSH). This can be seen by solving the model steady-state with carrying capacity only for the pituitary, = = 2 = 30 = = 0 , ≠ 0: Note that depends on many of the model parameters-= 3 3 , so that x2 depends now not only on thyrotroph production and removal rates, but also on TH production and removal rates, thyrocyte production and removal rates, and the carrying capacity of the thyroid gland.
Adding a carrying capacity term for the thyroid gland breaks the adaptation of x3 (TH). This can be seen by solving the model steady-state with carrying capacity only for the thyroid, = 30 = 2 = = = 0, ≠ 0: Again, KP depends on the model parameters, = 1 2 2 1 2 2 , making x3 sensitive to change in these parameters.

Computation of nullclines
To compute the nullclines / = 0, / = 0 we use separation of time scales, using the fact that hormone turnover times are much faster than the gland turnover times. We separate the model to the hormone "fast" equations (6)-(8), and the gland "slow" equations (9)-(10).
We first solve the steady-states for the hormone equations. The solution as a function of P and T can be represented as the roots of a third degree polynomial. For the sake of clarity, we express x1 and x3 using x2, and write an implicit equation for x2: With this, the "slow time scale" system is composed of two ODE's for P and T and one implicit algebraic equation for x2: In order to find the nullclines we need to solve the equations dT/dt=0 and dP/dt=0. One trivial solution is T=0, P=0.
In order to find the non-trivial solution we solve the equations T'=0 and P'=0 for x2, and substitute the solution into the implicit equation for x2. This gives the following nullclines: T'=0: The nullclines in the simple case (Fig. EV3A), i.e. 30 = 2 = = 0, takes the simple form of:

Stability analysis
In order to validate the stability of the slow system we need to find the signs of the eigenvalues of the Jacobian of the reduced system (P,T) (equations (11)-(13)) at the steady state.
We illustrate the different stability regimes in Figure EV3.
(i) We first inspect the solution at T>0, P>0 (Fig EV3B,C). To calculate the stability of this fixed point we need to calculate the partial derivative of x2 with respect to P and T.

Therefore, when this fixed point exists it is stable.
(ii) We next inspect the stability of the second fixed point at P=0, T>0 (Fig EV3D,E). This fixed point can be computed explicitly: In the simple case when there is no external thyroid hormone supply B30=0, and KX2=0, this condition is reduced to AB>1+KT.
(iii) Third, we inspect the fixed point at (P=0, T=0) ( Fig EV3F). This fixed point can be explicitly computed: The eigenvalues are negative providing that AB<1 and B30>1.
(iv) Finally, we inspect the stability of the fourth fixed point at T=0, P>0 (Fig EV3G). This point can be computed explicitly as:

Dependence on antibody parameter in Graves' disease
In this section we analytically explore the model dynamics under perturbation of the Ab parameter in Graves' diseases. We use the results of the stability analysis done above. In Graves' disease, the pituitary gland is far from its carrying capacity ≪ 1 , and there is no external thyroid hormone supply 30 = 0. We also assume for simplicity that 2 = 0. Therefore, the system has two fixed points in the P-T plane (Fig EV4A,B, Fig 4B): The parameter determines the sign and the stability of these two fixed points. When < 1, we get that only fixed point (i) is positive and stable (black full line in Fig EV4A,B). Above this value, when 1 < < 1 + , another unstable positive fixed point arises -fixed point (ii) (blue dashed line in Fig EV4A,B). Up to this point, the thyroid gland T has a constant value: = 1 1+ , which does not depend on . The value of thyroid hormones TH (x3) does not depend on as well (Fig EV4C). This is possible since the values of P and TSH (x2) do depend on AB -when AB is increased, P shrinks and TSH levels drop (Fig EV4D), compensating for the stimulatory effect of AB. However, when > 1 + , fixed point (i) becomes negative and unstable (black dashed line in Fig EV4A,B) while fixed point (ii) becomes stable (blue full line in Fig EV4A,B), a situation called a transcritical bifurcation. At this point, the pituitary gland has shrunk to zero, TSH=0, and the stimulatory effect of autoantibodies Ab cannot be compensated anymore. From this point on, thyroid gland size T and thyroid hormone level TH rise gradually with Ab: = 1 (1 − 1 ), 3 = −1 . These dynamics explain the occurrence of subclinical hyperthyroidism, in which TSH levels are low though thyroid hormone levels are kept normal (Fig 4B).

The relation between the steady-states of TSH and T4
In this section we compute the relation between TSH and T4, here denoted x2 and x3, that appears in Figure 4. Note that in this model, we ignore the complexity of T4 and T3 relation, and treat them as one variable x3.
From the steady-state of the equations for x1, x2 and P (1),(2),(5) we get: Note that 2 0 , 3 0 are not the steady-state for the full system with carrying capacities.
Assuming that in the normal healthy set-point the glands are far from their carrying capacities, we can thus compute , . We consider a normal set-point of FT4 3 0 = 15 / , and of TSH 2 0 = 1.5 / , and we normalize the pituitary mass units so that in the normal set point = 0 5 following measurements from Khawaja et al. that showed reduction of up to 80% in pituitary volume in hypothyroid patients following treatment. We thus derived a nonparameterized curve for x2(x3). Comparing this curve with data from Midgley et al (Midgley et al., 2013) yield a reasonable agreement (Fig 4D).

Changes in fetal thyroid and pituitary volume during gestation
A recent study found both pituitary and thyroid volume to correlate positively to gestational age in preterm and term infants but to correlate negatively to chronological age (Otani et al., 2021). Changes in volume during gestation may result from many different physiological effects. If we assume that the main change is due to the thyroid axis modeled here, the coordinated changes in P and T offer interesting constraints on which parameter group changes during gestation.
The thyroid and pituitary set points in the model, assuming they are not close to their carrying capacities, are: The change in T,P set-point before and after birth suggested by Otani et al. can thus reflect change in parameters in utero and after birth. Interestingly, P and T are observed to change together: they grow together before birth, and shrink together after birth. This fact sets a constraint on which parameter is responsible for this transition, because this parameter should affect P and T in the same direction. Therefore, the model combined with the findings of Otani et al. suggests that the parameter responsible for the opposite trends in gland volume change before and after birth is either and/or : the pituitary thyrotroph production and removal rates. This would be the case if thyrotroph proliferation rate is higher in utero than after birth. Note that these pituitary parameters dictate the volume accumulation of both glands -pituitary and thyroid. Another interesting prediction is the pituitary volume is more sensitive to changes in and due to the quadratic dependence on and , compared to the linear dependency of T on these parameters. This is in line with Otani et al that showed larger relative volume changes for P than for T for the same time period.

Trans-differentiation from other cell types into thyrotrophs
In this section, we consider other mechanisms for thyrotroph mass accumulation other than thyrotroph proliferation or growth. This includes trans-differentiation from other pituitary cell types (Horvath et al., 1990;Radian et al., 2003) and/or redifferentiation from dedifferentiated thyrotrophs (Wang et al., 2014).
Mathematically, we account for these effects by adding to the equation that describes the pituitary thyrotroph mass dynamics a source term inhibited by thyroid hormone levels. (16) The other model equations (1)-(4) remain the same.
We find that adding this term leaves the results of this analysis qualitatively similar.
As before, we rescale the variables by their steady-state value in the simple case in which the glands are far from their carrying capacities (i.e. infinite carrying capacity kT=kP=0), there are no TSHR antibodies and no external thyroid hormone supply (Ab=0, b30=0), and x3 depends approximately linearly on x2 (kx2=0). However here we keep the parameter that describes the thyrotroph proliferation/growth rate in the equations, so latter it can be compared to , the rate of thyrotroph mass accumulation from non-thyrotroph cells: The new scaled equations are:  Computing the steady-state for this system in the simple case described above, we find that the steady-state of x3 (thyroid hormone) is now dependent on the thyrotroph mass growth rate , and on the trans-differentiation term ( Figure S1A) . However, if the sum of rates + remains constant, x3 does not change. In other words, the overall rate of thyrotroph mass accumulation matters, but not their specific source.
When P approaches its carrying capacity, i.e. > 0, the symmetry between the transdifferentiation model to the original model starts to break ( Figure S1B). In this case, replacing thyrotroph proliferation bp with thyrotroph synthesis from non-thyrotroph source , has a similar effect on thyroid hormone levels but with scaling: Higher values map to lower bp values. Therefore, higher trans-differentiation rate is required to maintain a similar thyrotroph mass.
Next, we tested the effect of adding a trans-differentiation term on the TSH(TH) relationship ( Figure S1C). The models differ in their prediction for the TSH(TH) relation in hyperthyroidism. The original model predicts a sharp drop to very low TSH values even in subclinical hyperthyroidism with thyroid hormone levels that are only slightly above normal. However, the new model predicts a gradual drop in TSH with the rise in TH levels.
This difference stems from the behaviors of the two models in the regime of small thyrotroph mass. In the original model, when thyrotroph mass is equal to zero it cannot be restored again. Therefore, in Graves' disease, for high enough antibodies level, the model has a fixed point at zero thyrotroph mass: { 1 = −1 , 2 = 0, 3 = −1 , = 0, = 1 (1 − 1 )} (See Appendix section "Dependence on antibody parameter in Graves' disease" and Figure EV4). In this regime, the thyrotroph mass and TSH levels cannot compensate for the antibody stimulation, leading to increased TH levels and enlarged of thyroid mass. However, in the trans-differentiation model, thyrotrophs mass is generated even when it is low, due to the non-thyrotroph sources. Therefore, the fixed point at P=0 is lost. More data is needed to conclusively differentiate between these predictions.
Nullcline analysis and stream plot simulations reveal that the full dynamics, and not only the steady-state, remain qualitatively similar (Fig S1D-F). The nullclines for the model when 2 = 30 = = 0 are: However, a major difference between the models is that the original model has a nullcline at P=0, while the new model doesn't. This leads to the difference described above, in the regime of high thyroid mass T and small thyrotroph mass P.

Rise of TSH mRNA over days following thyroidectomy
Nolan et al provide evidence of a rise in TSH mRNA over days after thyroidectomy (Nolan et al., 2004). The present model captures this effect, based on the lifetime of TH. This lifetime is seven days in humans and is thought to be on the order of a few days in mice. Thus, after thyroidectomy, thyroid hormone levels fall gradually, reaching halfway down after its half-life. In this timescale of days, TSH mRNA expression in the thyrotrophs will be gradually de-repressed, and their levels would rise (Fig EV2). In contrast, thyrotroph mass will grow on a longer time scale of about a month. Figure EV2 illustrates this point: After thyroidectomy at time 0, T4 drops and TSH increases, reaching their minimal/maximal concentrations respectively after about 7 days. However, thyrotroph mass P does not reach its maximal point even after 30 days.